home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / RAND250.C < prev    next >
C/C++ Source or Header  |  1996-03-26  |  8KB  |  224 lines
  1. /*** Function prototypes for R250 random number generator, by W. L. Maier
  2. ***/
  3.  
  4. /* ------------------- */
  5. /* FUNCTION PROTOTYPES */
  6. /* ------------------- */
  7. # undef F
  8. # if defined(__STDC__) || defined(__PROTO__)
  9. #    define  F( P )  P
  10. # else
  11. #    define  F( P )  ()
  12. # endif
  13.  
  14. /* INDENT OFF */
  15. extern    void    Rand250_init F((int));
  16. extern    double    dRand250 F((void));
  17. extern    unsigned int Rand250 F((void));
  18. extern    unsigned int Rand250n F((unsigned int));
  19. static    unsigned int myrand F((void));
  20. static    void    mysrand F((unsigned int));
  21.  
  22. # undef F
  23. /* INDENT ON */
  24.  
  25.  
  26. /******************************************************************************
  27. *  Module:  Rand250.c
  28. *  Description: implements R250 random number generator, from S.
  29. *  Kirkpatrick and E.  Stoll, Journal of Computational Physics, 40, p.
  30. *  517 (1981).
  31. *  Written by:      W. L. Maier
  32. *
  33. *    METHOD...
  34. *        16 parallel copies of a linear shift register with
  35. *        period 2^250 - 1.  FAR longer period than the usual
  36. *        linear congruent generator, and commonly faster as
  37. *        well.  (For details see the above paper, and the
  38. *        article in DDJ referenced below.)
  39. *
  40. *    HISTORY...
  41. *        Sep 92: Number returned by dRand250() is in range [0,1) instead
  42. *            of [0,1], so for example a random angle in the
  43. *            interval [0, 2*PI) can be calculated
  44. *            conveniently.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
  45. *        Aug 92: Initialization is optional.  Default condition is
  46. *            equivalent to initializing with the seed 12345,
  47. *            so that the sequence of random numbers begins:
  48. *            1173, 53403, 52352, 35341...  (9996 numbers
  49. *            skipped) ...57769, 14511, 46930, 11942, 7978,
  50. *            56163, 46506, 45768, 21162, 43113...  Using ^=
  51. *            operator rather than two separate statements.
  52. *            Initializing with own linear congruent
  53. *            pseudorandom number generator for portability.
  54. *            Function prototypes moved to a header file.
  55. *            Implemented Rand250n() to generate numbers
  56. *            uniformly distributed in a specific range
  57. *            [0,n), because the common expedient rand()%n is
  58. *            incorrect.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
  59. *        May 91: Published by W. L. Maier, "A Fast Pseudo Random Number
  60. *            Generator", Dr. Dobb's Journal #176.
  61.  
  62. ******************************************************************************/
  63.  
  64. #include <stdlib.h>
  65.  
  66. #define    NEXT_RAND(r250_index, K, L, M, new_rand)        \
  67. {                                \
  68.     register int j;                        \
  69.     if (r250_index >= K)                    \
  70.     j = r250_index - K;    /* Wrap pointer around */    \
  71.     else                            \
  72.     j = r250_index + L;                    \
  73.     new_rand = r250_buffer[r250_index] ^= r250_buffer[j];    \
  74.     if (r250_index >= M)    /* Increment for next time */    \
  75.     r250_index = 0;                                         \
  76.     else                                                        \
  77.     r250_index++;                                           \
  78. }
  79. /**** Static variables ****/
  80. static int r250_index = 0;
  81. static unsigned int r250_buffer[250] =
  82. {
  83.     15301, 57764, 10921, 56345, 19316, 43154, 54727, 49252, 32360, 49582,
  84.     26124, 25833, 34404, 11030, 26232, 13965, 16051, 63635, 55860, 5184,
  85.     15931, 39782, 16845, 11371, 38624, 10328, 9139, 1684, 48668, 59388,
  86.     13297, 1364, 56028, 15687, 63279, 27771, 5277, 44628, 31973, 46977,
  87.     16327, 23408, 36065, 52272, 33610, 61549, 58364, 3472, 21367, 56357,
  88.     56345, 54035, 7712, 55884, 39774, 10241, 50164, 47995, 1718, 46887,
  89.     47892, 6010, 29575, 54972, 30458, 21966, 54449, 10387, 4492, 644,
  90.     57031, 41607, 61820, 54588, 40849, 54052, 59875, 43128, 50370, 44691,
  91.     286, 12071, 3574, 61384, 15592, 45677, 9711, 23022, 35256, 45493,
  92.     48913, 146, 9053, 5881, 36635, 43280, 53464, 8529, 34344, 64955,
  93.     38266, 12730, 101, 16208, 12607, 58921, 22036, 8221, 31337, 11984,
  94.     20290, 26734, 19552, 48, 31940, 43448, 34762, 53344, 60664, 12809,
  95.     57318, 17436, 44730, 19375, 30, 17425, 14117, 5416, 23853, 55783,
  96.     57995, 32074, 26526, 2192, 11447, 11, 53446, 35152, 64610, 64883,
  97.     26899, 25357, 7667, 3577, 39414, 51161, 4, 58427, 57342, 58557,
  98.     53233, 1066, 29237, 36808, 19370, 17493, 37568, 3, 61468, 38876,
  99.     17586, 64937, 21716, 56472, 58160, 44955, 55221, 63880, 1, 32200,
  100.     62066, 22911, 24090, 10438, 40783, 36364, 14999, 2489, 43284, 9898,
  101.     39612, 9245, 593, 34857, 41054, 30162, 65497, 53340, 27209, 45417,
  102.     37497, 4612, 58397, 52910, 56313, 62716, 22377, 40310, 15190, 34471,
  103.     64005, 18090, 11326, 50839, 62901, 59284, 5580, 15231, 9467, 13161,
  104.     58500, 7259, 317, 50968, 2962, 23006, 32280, 6994, 18751, 5148,
  105.     52739, 49370, 51892, 18552, 52264, 54031, 2804, 17360, 1919, 19639,
  106.     2323, 9448, 43821, 11022, 45500, 31509, 49180, 35598, 38883, 19754,
  107.     987, 11521, 55494, 38056, 20664, 2629, 50986, 31009, 54043, 59743
  108. };
  109.  
  110. /**** Function: Rand250_init
  111.     Description: initializes Rand250 random number generator. ****/
  112.  
  113. void
  114. Rand250_init(int seed)
  115. {
  116.  
  117. /*---------------------------------------------------------------------------*/
  118.     int     j, k;
  119.     unsigned int mask;
  120.     unsigned int msb;
  121.  
  122. /*---------------------------------------------------------------------------*/
  123.     mysrand(seed);
  124.     r250_index = 0;
  125.     /* Fill the r250 buffer with 15-bit values */
  126.     for (j = 0; j < 250; j++)
  127.     r250_buffer[j] = myrand();
  128.     for (j = 0; j < 250; j++)        /* Set some of the MS bits to 1    */
  129.     if (myrand() > 16384)
  130.         r250_buffer[j] |= 0x8000;
  131.     msb = 0x8000;            /* To turn on the diagonal bit    */
  132.     mask = 0xffff;            /* To turn off leftmost bits    */
  133.     for (j = 0; j < 16; j++)
  134.     {
  135.     k = 11 * j + 3;            /* Select a word to operate on  */
  136.     r250_buffer[k] &= mask;        /* Turn off bits left of diag.    */
  137.     r250_buffer[k] |= msb;        /* Turn on the diagonal bit     */
  138.     mask >>= 1;
  139.     msb >>= 1;
  140.     }
  141. }
  142.  
  143. /**** Function: Rand250 Description: returns a random unsigned integer k
  144.         uniformly distributed in the interval 0 <= k < 65536.  ****/
  145.  
  146. unsigned int
  147. Rand250(void)
  148. {
  149.  
  150. /*---------------------------------------------------------------------------*/
  151.     register unsigned int new_rand;
  152.  
  153. /*---------------------------------------------------------------------------*/
  154.  
  155.     NEXT_RAND(r250_index, 147, 103, 249, new_rand);
  156.  
  157.     return new_rand;
  158. }
  159.  
  160. /**** Function: Rand250n Description:
  161.  
  162. returns a random unsigned integer k uniformly distributed in the interval
  163.     0 <= k < n    ****/
  164.  
  165. unsigned int
  166. Rand250n(unsigned int n)
  167. {
  168.  
  169. /*---------------------------------------------------------------------------*/
  170.     register unsigned int new_rand, limit;
  171.  
  172. /*---------------------------------------------------------------------------*/
  173.     limit = (65535U / n) * n;
  174.  
  175.     do
  176.     {
  177.     NEXT_RAND(r250_index, 147, 103, 249, new_rand);
  178.     }
  179.     while (new_rand >= limit);
  180.  
  181.     return (new_rand % n);
  182. }
  183.  
  184. /**** Function: dRand250
  185.     Description: returns a random double z in range 0 <= z < 1.  ****/
  186. double
  187. dRand250(void)
  188. {
  189.  
  190. /*---------------------------------------------------------------------------*/
  191.     register unsigned int new_rand;
  192.  
  193. /*---------------------------------------------------------------------------*/
  194.  
  195.     NEXT_RAND(r250_index, 147, 103, 249, new_rand);
  196.  
  197.     return new_rand / 65536.;        /* Return a number in [0.0 to 1.0) */
  198. }
  199.  
  200.  
  201. /***  linear congruent pseudorandom number generator for initialization ***/
  202.  
  203. static unsigned long seed = 1;
  204.  
  205.  /*
  206.   * Return a pseudorandom number in the interval 0 <= n < 32768. This
  207.   * produces the following sequence of pseudorandom numbers: 346, 130,
  208.   * 10982, 1090...  (9996 numbers skipped) ...23369, 2020, 5703, 12762,
  209.   * 10828, 16252, 28648, 27041, 23444, 6604...
  210.   */
  211. static unsigned int
  212. myrand(void)
  213. {
  214.     seed = seed * 0x015a4e35L + 1;
  215.     return(unsigned) (seed >> 16) & 0x7fff;
  216. }
  217.  
  218. /*    Initialize above linear congruent pseudorandom number generator */
  219. static void
  220. mysrand(unsigned int newseed)
  221. {
  222.     seed = newseed;
  223. }
  224.